Numerical detection of stochastic to deterministic transition 



R.K. Brojen Singh 

Centre for Interdisciplinary Research m Basic Sciences, Jamia Millia Islamia,New Delhi 110025, India. 
* Corresponding Author at: Centre for Interdisciplinary Research in Basic Sciences, 
Jamia Millia Islamia, New Delhi 110025, India, 
Tel: +911-2698 1111 (4492) Telfax: +911-2698 3409 and 
Email:rksingh@]mi. ac. in 
(Dated: September 2, 2011) 

ABSTRACT 

We present the numerical estimation of noise parameter induced in the dynamics of the variables 
by random particle interactions involved in the stochastic chemical oscillator and use it as order 
parameter to detect the transition from stochastic to deterministic regime. In stochastic regime, 
this noise parameter is found to be increased as system size decreases, whereas, in deterministic 
regime it remains constant to minimum value as system size increases. This let the transition 
from fluctuating to fixed limit cycle oscillation as the system goes from stochastic to deterministic 
transition. We also numerically estimated the strength of the noise parameter involved both in 
Chemical Langevin equation and Master equation formalisms and found that strength of this 
parameter is much smaller in the former than the later. 

ifEYWORDS : Noise parameter. Chemical oscillator, fluctuating limit cycle oscillation, Master 
equation. Chemical Langevin equation. 
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The origin of noise (intrinsic and extrinsic) in stochas- 
tic systems is beheved to be due to random particle or 
molecular interactions taking place among the particles 
or molecules in the physical, chemical or biological sys- 
tems and various forms of fluctuations in the environment 
surrounding the system [l|, 0] . The strength of the noise 
in such complex systems depends on various parameters, 
for example, system size V, population of particles or 
molecules accomodated in the systems {N) and various 
types of fluctuations in the surrounding environment etc 
0, 01 ■ For instance, the noise strength associated with 
single cell gene expression scales as of relative fluc- 
tuation amplitude [sj . The estimation of the noise associ- 
ated with stochastic systems via Master equation formal- 
ism could tell us the role of noise in the system dynamics 
0, Q • For instance, noise become an essential parameter 
which plays a constructive role in biological systems and 
is used in weak signal amplification and enhancing the 
detection of information carrying weak signal, the phe- 
nomenon known as stochastic resonance 0, . Further 
this parameter might also be used by a group of systems 
as a means of synchronizing or correlating behaviors of 
each individual systems as is seen in various multicellu- 
lar organisms Q and in a group of identical unicellular 
organisms (Toj . 

The deterministic systems can be considered as noise 
free systems 0, [ll| . Because the noise associated with 
the dynamics of each variables involved in determinis- 
tic systems becomes negligible [ll|. The transition from 
stochastic to deterministic system can be scaled by defin- 
ing a thermodynamics limit, 1^ — >■ oo, iV — > oo but 
N/V —7- finite [1, Attempts have been made in 

genetic oscillator by calculating noise-to-signal ratio as a 



function of some reaction rate to differentiate determin- 
istic from stochastic systems [l^j- However estimation 
of this limit or to define correct order parameter to dis- 
tinguish between stochastic and deterministic systems is 
still an open question both for numerical as well as ex- 
perimental situations. 

The noise fiuctuation is an inherent property of most 
of the biological systems and has two contrast roles de- 
pending upon its magnitude. If the strength of the noise 
is large enough (above a critical limit, say ifjijc), then it 
destructs the signal. However if the strength of the noise 
is weak, say 77(770 then its role becomes constructive and 
is used to detect and amplify the weak signal (stochastic 
resonance) [1, Q . In this work we try to identify a pa- 
rameter that can separate stochastic and deterministic 
systems. Probably the study may unfold many inter- 
esting roles of noise in various biological systems. For 
this purpose we study a chemical oscillator defined by an 
oregonator reaction network via Master equation formal- 
ism. Then we briefiy explain the stochastic simulation 
algorithm (SSA) Q to be implemented to simulate the 
reaction network. Further we study Chemical Langevin 
formalism of this model to measure noise parameter in- 
duced in the system dynamics. It is followed by results 
and discussions, and then few conclusions based on the 
results we have got. 

The intrinsic noise in stochastic systems is due to ran- 
dom particle or molecular interactions taking place in the 
system and is an inherent property associated with the 
stochastic system Q. If we consider a configurational 
state x{t) = [xi,X2, ■■■,xn]~^ consisting of N species 
(molecule or particle) at any instant of time t, then the 
random interaction events (only collisions which gives rise 
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FIG. 1: The plots of the dynamics oi X, Y amd Z at three 
different system sizes, V = 1, 100 and 1000 respectively simu- 
lated using Gillespie's SSA [3]. The parameters used are same 
as IS used in Gillespie, (1977) [1]. 



decay and creation of particle) taking place in the system 
is given by, 



SqXg + ■ 



(1) 



where, {S} are co-efficients of the species in the /xth in- 
teraction. kfj_, with n = 1,2, ■ ■ ■ ,Mis the rate of interac- 
tion of the species {xi}, i — a,b, - ■ ■ to give species {xj}, 
j = p,q, ■ • ■ in a certain interval of time [t, t + A<] . This 
leads to the change in states from x{t) at time t to a new 
state x'{t+At) during the time interval. Depending upon 
the magnitude of At, there could be L number of interac- 
tions in series took placed. If we simplify the transitions 
by taking an infinitesimal time interval dt such that dur- 
ing [t, t -\- dt] only one interaction is occured then we can 
write At = dtl^l -|- dtl^l, • • • , +dt^^\ which leads to a se- 
ries of microscopic state changes, x x'^' — > .f 1^1 • • • — > 
a;[^] , x' . These {dt^^ , j = 1,2, L are not necessarily 
the same but taking dtl^l = d^I^l = • • • d^I^l = dt, we have 

At ~ Ldt leading to macroscopic state change x x' . 
Now if we consider each microscopic state change dur- 



ing [t + dtbl , t + , {j, j' = 1, 2, • • • , L; j ^ /} then 
the time evolution of the state probability, Pj{xj,t) of 
the state change based on transition probabilitis {W} of 
decay and creation of the particles evolved in the inter- 
action event can be described by the following Master 
equation 0,01, 



9P3iXj,t) 

dt 



P,{x„t)Ws,- 



^Y,P,{x-,.,t)Ws^, 



(2) 



The macroscopic state change can be described by a se- 
ries of state probabilities corresponding to each micro- 
scopic state changes P{x, t) Pi ^ P2 ■ ■ ■ Pl ^ 
P{x',t). Since each interaction that drives a particular 
microscopic state change is random in nature, the tra- 
jectory of Pj from P{x,t) to P{x',t) will a Brownian 
trajectory. 

We consider the chemical oscillator model known as 
orcgonator devised by Field and Noves |l3| based on the 
criticism made by Tyson and Light [ij] on Brusselator 
model which is two molecular species reaction model. 
The modified chemical oscillator model consists of three 
molecular species, X, Y and Z involved in the following 
five reaction channels. 



A + Y ^ X 

X + Y H B 

C + X H 2X + Z 

k4 



(3) 
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where {ki}, i = 1,2,-- - ,5 are reaction rate constants 
and A, B, C, D and E are constants. The state of 
the system at any instant of time t is given by the vec- 
tor, x{t) = [X{t),Y{t), Z{t)]~^. Considering only micro- 
scopic state change, we can construct molecular trans- 
formation diagram [7| and derive transition probabilities 
{W} of each reaction channel during the time interval 
[t,t + dt]. Depending upon the arrow of the diagram 
and state change due to molecular decays and creations 
involved in the molecular interactions indicated by the 
reactions listed in|31 we obtain the following Master equa- 
tion (ME) of the system, 



-PiX,Y,Z;t) 



= kiA{Y + 1)P{X -l,Y + l,Z;t) + k2{X + 1){Y + 1)P(X + l,Y + l,Z;t) 



+k3C{X - 1)P{X -1,Y,Z-1) + -kiX{X + 1)P{X + 1, Y, Z) + ksEiZ + 1)P{X, Y-1,Z + 1) 

(4) 



[kiAY + k2XY + k^CX + kiX"^ + k^EZ] P{X, Y, Z) 
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FIG. 2: The plots showing the estimation of noise parameter in amphtudes for X and Y variables. The left panels show the 
transition from fluctuating to fixed limit cycle oscillation as the function of V, as V goes from small (stochastic regime) to 
large V (deterministic regime). The right panels show the phase plot in {tja, V) plane indicating stochastic and deterministic 
regimes. 



However, it is very difficult to solve this Master equa- 
tion (U). The alternative way to solve such complicated 
molecular processes in stochastic system given by this 
equation is to do simulation by simplifying the process of 
jumping from one stochastic state to another as discrete 
Markov process Q . This stochastic simulation algorithm 
(SSA) due to Gillespie, is based on the assumption that 
the time interval to jump from one stochastic state to an- 
other is taken to be small enough such that at most one 
reaction type can occur. The algorithm systematically 
allows to pick up each reaction event randomly from the 
reaction set at every time step to allow transition from 
one state to another along the trajectory of the of the 
variables by defining a joint probability density function, 
r(r, /Lt) = n(T)A(/z). A(/z) is the probability density func- 
tion of picking up /xth reaction and n(r) is the probabil- 
ity density function that at time step r the reaction ji 
will fire. The reaction time and reaction number fired at 
that time can be estimated computationally by generat- 
ing two uniform random numbers ri and r2 to identify 



T and by T = ^-ZJ-^'^^ 77 and = r2j2 

tively by imposing the condition J2z=i ^ fuiYmi 
where {w} are propensity functions. 
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Now we solve the stochastic system described by the 
orcgonator reaction model given by ^ using SSA and 
the results are shown in Fig.l. In the figure we present 
the dynamics of X, Y and Z as a function of time for 
different values of system sizes V. It is seen from the dy- 
namics that at small V{=1), the noise fluctuation associ- 
ated with each variable dynamics is large as compared to 
the behavior at large F(=100, 1000). Further the ampli- 
tudes of each variable oscillation at small V are random 
in nature due to noise. The degree of randomness in am- 
plitudes which is proportional to the magnitude of noise 
is indicated by the width of the fluctuating limit cycle 
oscillation shown in Fig. 2 left panels (plots in XZ and 
YZ planes). However as V increases the amplitudes be- 
come ordered and move towards a fixed value tending to 
limit cycle oscillation at that fixed value of amplitude as 
is shown in Fig. 2. This disorder to orderness in ampli- 
tudes as a function of V could be one way to look at 
the transition from stochastic to deterministic system by 
defining an order parameter 77^ which is noise in ampli- 
tude given by, tja = j^- Here {A) is the mean of the 
random amplitudes of the spikes for a particular V and 
GA = (((^) — AY)^/'^ is the standard deviation of the 
amplitudes. The transition from stochastic to determin- 




istic can be identified by looking at the behavior of ry^ as 
a function of V. At thermodynamic hmit i.e. at ^ — S- oo, 
essentially rjA because of the transition from disor- 
der or randomness to orderness in A. Numerically one 
can define a critical value Vc such that: (1) if V{Vc, rjA 
will have values larger than a, -q^ a.sV decreases and so 
the system is in stochastic regime, (2) if V)Vc the values 
of rjA will remain almost stationary to rj^ as V increases, 
and therefore the system is in deterministic regime. 

We now simulate the oregonator model ^ to obtain 
r]A for each value of V by averaging the amplitudes of 
300 spikes for each value of V . The plots of rjA for X 
and Y variables are shown in the right hand panels in 
Fig. 2. From these plots we could able to identify the 
approximate critical values of V and rjA to heVc^ 70 ± 5 
and Tjc ^ 0.004 ± 0.001 respectively. In the stochastic 
regime {V{Vc) tja increases as V decreases, whereas in 



V)Vc, r]A is stationary as V increases. The stochastic and 
deterministic regimes are shown in the {rjA — V) phase 
diagrams as shown in the panels. 

We now follow Gillespie's technique Q to reduce the 
Master equation (j4]) to a more simpler form, the Chemical 
Langevin Equation (CLE), based on two important real- 
istic approximations made on a random variable K{x, r) 
which is the number of a particular reaction fired during 
a time interval [t, t + r] with t)0. The first approximation 
is to impose small r limit that let the propensity function 
remain constant (ui ~ constant) during the time interval 
and K to approximate to statistically independent Pois- 
son distribution functions. The second approximation is 
to impose large r limit that let and Poisson distri- 
bution function to be replaced by Normal distribution 
function with same mean and variance. Following these 
steps we reach the following CLE for oregonator. 



dX 
'dt 
dY 
~dt 
dZ 
H 



kiAY - k2XY + ksCX - k^X^ 

-kiAY - k2XY + k^EZ ^ 
1 



VkiAY^i - Vk2XY^2 + VhCX^s - ^kJC^^i 



WkiAYi^ - Vk2XY^e + ^k^EZ^y 



(5) 



fcgCX - k^EZ 



VksCX^s - VhEZ^9 
I 



where, = lim.dt^()Ni{Q,l) / \/dt are noise parameters 
satisfying, £,i{t)^j{t') = S,jS{t - t'). {A^(0, 1)} are nor- 
mal distribution functions with mean and variance 1. 
The noise terms in the dynamics of each variables are 



functions of V, £, and the variables in the system. If the 
noise term is negligible, the equation ^ reduces to de- 
terministic equations whose steady state solutions can be 
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obtained by putting -^x{t) = which are given by, 
_ hA ( I 8Ck,k, \ 

4fc| 1^ Akik^ V Akiki J ^ ' 

^ hksAcf j 8Ck2k3 \ 
2k2kzE yl Akiki j 

If wc represent Xc{t) and Xd{t) as the variables obtained 
by solving CLE in ([S]) and deterministic equations, then 
the noise term associated with variables Xc{i) can be es- 
timated by rjL[t) ^ [xc{t) -~Xd[t)\. From equation (O we 
can write the noise term as, 

77L(fd, lV;t) ^ j [xc{t) - Xd{t)] dt 

= ^ j F{xd,lV-t)dt (7) 

where F{xdT£,,V]t) is a function and 

^=[^i,^2,''- iCw]"^! iV = 9 in the case of oregona- 
tor model we consider. This noise parameter can also be 
used to detect stochastic to deterministic transition in 
the same fashion we have done in the case of stochastic 
simulation. Now we solve CLE ([5]) using standard 
Runge-Kutta method for numerical integration of a set 
of differential equations [l^ for different values of V . 
For each V we calculate noise parameter in amplitude 
•qL as we have done above in the case of SSA by taking 
300 spikes. The results are shown in Fig. 3. The results 
show that r]^{rj^))rjA{X,Y) in magnitudes. This shows 
that CLE system described by ([5]) is much much closer 
to deterministic system as compared to ME system 
given by (g]). 

Since there is no exact and clear demarcation between 
actual stochastic and deterministic systems, the attempt 
to detect the transition between these two regimes with 
exact number is almost not possible. However we can 
look for some order parameter which can able to detect 
this transition approximately. For example we used noise 
parameter, which is calculated from random amplitudes 
as a function of system size, as order parameter to de- 
tect this transition. The randomness in the amplitudes 
of each individual spikes in the dynamics of the variables 
involved in a stochastic system in a particular system size 
is due to noise associated with these variables induced by 



random particle interaction in the system. This random- 
ness in amplitudes will become order and tends to a fixed 
value of amplitude as the stochastic system iVs) goes to 
deterministic limit (V^) i.e as [Vs Vd)- This leads to 
the transition from the fluctuating limit cycle oscillations 
to the fixed limit cycle oscillation. This transition can be 
well characterized by the noise parameter in amplitudes 
{va, Vl) and we use it as an order parameter to detect 
the transition. We could able to detect the transition 
from the behavior of r]A which shows increasing nature 
as V decreases and minimum stationary value as V in- 
creases. However there are some issues in our numerical 
experiments, for example it is very difficult to find ideal 
deterministic limit where rjA should be zero. Further it 
is hard to decide the number of spikes to be taken for 
calculating (A) (we took 300 spikes) so that the value of 
r]A to be obtained is statistically significant and correct. 

The Master equation formalism of the chemical oscil- 
lator model we have taken and the SSA simulation sys- 
tematically takes into account the noise fluctuations in 
the dynamics of the variables from the particle interac- 
tion picture. However, the Chemical Langevin Equation 
beautifully can able to connect stochastic and determin- 
istic descriptions scaling by volume and the noise terms 
fluctuating with the order of 1/"^/^ letting the variables 
evolve with time. But the value of r)^ and rj^ obtained 
in CLE is found to be much smaller than the value rjA 
found in ME. 

Most living systems probably use noise fluctuations 
in a more constructive way, within an optitimal level, 
in weak signal amplification and information processing. 
However if the strength of the noise is larger than this 
optimal value then the role of the noise become in the 
destruction of the signal. It is also to be noted that the 
size of living cells from birth to death also fluctuate in 
short time scales. In such situations the role of noise 
could be different from time to time whenever changes in 
the cellular size takes place. Such changes in the role of 
noise could affect on various biological functions, cellular 
communications, signal manipulations etc. The investi- 
gation of such problems could explore important issues 
regarding how biological, chemical and physical systems 
work. 
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